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A SABR-CODED FAST FOURIER TRANSFORM SUBROUTINE 


Gerald N. Cederquist 

Cooley Electronics Laboratory 
University of Michigan 
Ann Arbor, Michigan 48105 

ABSTRACT 

Using the Fast Fourier Transform algorithm, this subrou¬ 
tine computes in situ either the direct or inverse discrete Fourier 
transform of a pure-power-of-two number of complex points. 
Written for use with the DEC 8K FORTRAN system, it runs in less 
time and takes less core space than the same algorithm coded in 
FORTRAN. It has been extensively tested and checked; for example, 
the subroutine will do a direct followed by an inverse transform of 
128 points in 47 seconds with a round trip root mean square error 
of 1.23 parts per million. 

INTRODUCTION 

The Acoustical Signal Processing Group at Cooley Labora¬ 
tory is engaged in data acquisition and analysis of underwater sound 
signals. Due to constantly changing requirements in our data analy¬ 
ses, we write most of our analysis programs in FORTRAN and use 










the DEC 8K FORTRAN System to run these programs. We have 
found the Fast Fourier Transform subroutine documented here to 
be a useful tool in data analysis and wish to pass it along to others 
who can make use of it. 

Note that by using FORTRAN, we are trading CPU time for 
the capability to make quick and easy modifications to our programs. 
We do not recommend using this FFT subroutine (or 8K FORTRAN 
for that matter) in computation-intensive programs; these programs 
should be written in assembly code or (preferably) run on a larger 
machine than the PDP- 8. 

REQUIREMENTS 

Hardware 

Required hardware to use the FFT subroutine is the same as 
the required hardware to use the DEC 8K FORTRAN system: a pro¬ 
cessor which will execute the PDP-8 order code and which has at 
least 8K of memory. Thus an 8K PDP-8 or an 8K LINC-8 or an 
8K PDP-12 may be used. Note that the subroutine does no I/O and 
thus does not place any requirement on the I/O devices attached to 
the machine. 

Storage 


The FFT subroutine requires 5 pages of memory in any 









memory bank to store the instructions for computing the FFT. In 
addition the subroutine requires autoindex register 10 and locations 
100-120 inclusive on page 0 of the memory bank in which FFT re¬ 
sides . 

Note that the subroutine does not itself contain the data area 
holding the points to be transformed. Thus it is possible to trans¬ 
form up to 512 points (stored in COMMON) in an 8K machine and 1024 
points in a 12K or larger machine. 

Subroutines 

The following subroutines from the 8K FORTRAN Library 
are required: 

SIN, COS, FAD, FSB, FMP, FDV, FLOT, STO, CHS 

USAGE 

Loading 

The FFT subroutine must be assembled using the SABR 
assembler. (See DEC manual DEC-08-ARXA-D, 8K SABR Assembler , 
available from the DEC Program Library.) The resulting binary 
tape may be loaded by the 8K Relocating Loader into any program 
which needs the FFT subroutine. 
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Calling Sequence 

The FFT subroutine may be invoked via the FORTRAN state¬ 


ment 


CALL FFT (A, M, ISIGN) 


to compute either the direct transform when ISIGN is < 0: 


F, = 


2 M 

-4j-= Yj A. exp(-27ri • jk/2 M ), 
2 M j=l 1 


k = 1, ..2 


M 


or the inverse transform when ISIGN is > 0: 


>M 


JVL 


F = Yj a - ex P (2^i • jk/2 ) : 
k j=i J 


k = 1, ..2 


M 


where i = • (Though not indicated above, the transform is done 

in-situ, and needs no scratch storage for F .) 

The calling parameters have the following significance: 

"A” is an array of complex numbers to be transformed, 
stored in the format ’’real part, imaginary part, 
real part, imaginary part, ... The numbers 
forming each complex pair are assumed to be in 
FORTRAN floating point form. 

"M" is the logarithm to the base 2 of the number of 

points to be transformed, and must be a FORTRAN 




integer. 
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"ISIGN" is any positive integer if the sign on the exponent 
of the complex roots of unity to be used in the 
transform is to be positive, and is any negative 
integer if the sign on the exponent of the complex 
roots of unity is to be negative. Note that if ISIGN 
is negative, the transform will be taken, and then 
additionally, all values in the A array will be di¬ 
vided by 2**M (normalizing the direct transform). 

For example, the direct transform of a 256 point complex¬ 
valued time series stored in the array "X" may be taken by the 
statement 

CALL FFT(X, 8,-1) 

and the inverse transform of a 32 point complex spectrum stored in 
the array "SPECT" may be taken by the statement 

CALL FFT(SPECT, 5,1) 


EXECUTION TIME 

The execution time of the FFT subroutine for both direct and 
inverse transforms is given in the table below. Note that the fourth 
column contains an estimate of the time needed to transform accord¬ 
ing to the definition of the Discrete Fourier Transform. The 


I 
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difference is the time saved by using the FFT technique. 


Number of 
Points 

Direct 

FFT 

Inverse 

FFT 

By The 
Definition 

4 

< 1 

< 1 

1 

8 

1 

1 

3 

16 

2 

2 

8 

32 

5 

4 

29 

64 

11 

10 

112 (1.8 min) 

128 

24 

23 

430 (7.2 min) 

256 

53 

51 

1664 (27.7 min) 

512 

118 

113 

6571 (109.5 min) 


Table 1. Execution Time of FFT in seconds 

Note that it is only the FFT algorithm which makes computing dis¬ 
crete Fourier transforms feasible in 8K FORTRAN. 

ERROR ANALYSIS 

The FFT subroutine was extensively tested for the 128 point 
case on a time series of samples from a white Gaussian random pro¬ 
cess. The method used was to transform a time series into the 
frequency domain and back to time domain using FFT, and then to 
compute the root mean square error incurred during the round trip. 
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If the original time series is denoted by A i and the doubly-trans¬ 
formed time series by B., then the RMS error E is defined to be 


f 128 V 

L ( A i - B i) 2 


E = 


i=l 


V 


128 

2 V 


i=l 


/ 


as suggested by the discussion in Ref. 2. 

The value of E averaged over 512 trials on an integer valued 
time series having standard deviation of 249.6 in both real and 
imaginary components was 1.2380 x 10 . Increasing the standard 

deviation to slightly over one million leads to an E value of 1.2311 x 

g 

10 for 256 trials. In either test, the maximum value of E did not 
__ 0 

exceed 1.47 x 10 in any one trial. 

The value of the direct transform via the FFT subroutine was 
compared with the value obtained by analytical computation for the 
case of the 128 point waveform shown below: 




|.0 


0. 




1 


/ Ze.ro 


6* 


12 ? 



2.CTO gi/ery cohere. 


64 


IZf 
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_7 

The largest error found was 1.4x10 

As a final test, the average D of the difference between a 
set of original and doubly-transformed points was computed for 
M = 2 through 9: 




D = 


2 M 

2 

i=l 


(A; 






The table below gives the values of D obtained on the previously- 
mentioned white Gaussian process with standard deviation 249.6: 


M 



D 


2 

3 

4 

5 

6 

7 

8 
9 


4 

8 

16 

32 

64 

128 

256 

512 


0.E0 

1.90735E-6 
-2.86102E-6 
-2.38418E-7 
-4.76836E-7 
1.90734E-6 
-2.38418E-6 

l.E-5 (Upper bound. True value 
not known due to memory 
limitations) 


The conclusion we draw from all these tests is that the FFT 
subroutine is capable of giving slightly less than 6 figure accuracy 
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over the range M = 2 through M = 8 and at least 5 figure accuracy 
for M = 9 . 

LISTING 

The listing which follows this writeup was done on version 16 
of the SABR assembler running within CPS, the local Cooley Labor¬ 
er atory operating system for the LINC-8 (see DECUS program L-80). 

Total assembly time excluding listing is 48 seconds. 


REFERENCES 


J. W. Cooley, P. A. Lewis and P. D. Welch, "The Fast 
Fourier Transform and Its Applications, ” IEEE Transac- 
tions on Education , Vol. 12, No. 1, pp. 27'ff, March 1969. 
The algorithm used in the SABR-Coded FFT described here 
is a direct copy of the one shown in FORTRAN IV on page 29 
of this paper. 

P. D. Welch, "A Fixed-Point Fast Fourier Transform Error 
Analysis, " IEEE Transactions on Audio, Vol. 17. No. 2. 
p. 151, June 1969. 














PAGE 01 SABR ASSEMBLED V(16> 


-10- 


/SABR COOED FAST FOUR I ER TRANSFORM 
/SUBROUTINE FOR A PURE-POWER-OF-2 
/NUMBER OF COMPLEX POINTS 
/ 

/THIS SUBROUTINE FOR USE IN 8K FORTRAN PROGRAMS 

/GERALD N. CEDERQUIST 

/COOLEY ELECTRONICS LABORATORY 

/NORTH CAMPUS 

/UNIVERSITY OF MICHIGAN 

/ANN ARBOR* MICHIGAN 48105 

/I SEPTEMBER* 1969 

/ 

ENTRY FFT 
/ 

/CALLED vJITH PARAMETERS A* M* AND I SIGN 
/WHERE: 

/ # 

/"A" IS AN ARRAY OF COMPLEX NUMBERS TO BE 
/TRANSFORMED* STORED IN THE FORMAT 
/’’REAL PART* IMAGINARY PART*..." 

/THE NUMBERS FORMING EACH COMPLEX PAIR 
/ARE ASSUMED TO BE IN FORTRAN FLOATING 
/POINT FORM. 

/ 

/"M" IS THE LOGARITHM TO THE BASE 2 
/OF THE NUMBER OF POINTS TO BE 
/TRANSFORMED* AND MUST BE A FORTRAN INTEGER. 

/ 

/"ISIGN" IS ANY POSITIVE INTEGER IF 
/THE SIGN ON THE EXPONENT OF THE COMPLEX ROOTS OF 
/UNITY TO BE USED IN THE TRANSFORM IS TO 
/BE POSITIVE* AND IS ANY NEGATIVE INTEGER 
/IF THE SIGN ON THE EXPONENT OF THE COMPLEX ROOTS 
/OF UNITY IS TO BE NEGATIVE. NOTE THAT IF ISIGN 
/IS NEGATIVE* THE TRANSFORM WILL BE TAKEN* AND THEN 
/ADDITIONALLY* ALL VALUES IN THE A ARRAY WILL 
/BE DIVIDED BY 2**M 
/ 

/REFERENCE: COOLEY* LEWIS* AND WELCH* "THE 

/ FAST FOURIER TRANSFORM AND ITS APPLICATION^"* 

/ IEEE TRANSACTIONS ON EDUCATION* VOL 12* #1* 

/ PAGE 27 FF* MARCH* 1969. 

/ 

/STORAGE REQUIRED: 

/ AUTOINDEX REGISTER 10 

/ LOCNS 100-120 INCLUSIVE ON PAGE 0 OF THE 

/ MEMORY BANK IN WHICH FFT RESIDES 

/ 5 PAGES ABOVE 200C0CTAL) for the PROGRAM PROPER 

/ 

/SUBROUTINES REQUIRED: 

/ SIN* COS* FAD, FSB* FMP * FDV* FLOT* STO* CHS 

/ 

/NOTE THAT THE COMMENTS CONCERNING THE CODE FOR 
/THE FFT REFLECT THE FACT THAT IT WAS LIFTED FROM 
/THE FORTRAN PROGRAM GIVEN IN THE REFERENCE. 

/ALSO* THE FOLLOWING NOTATIONAL CONVENTIONS 
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/ARE USED: 

/ "->" MEANS "POINTER TO" 

/ .REAL. AND .IMAG. ARE THE REAL- AND 

/ IMAGINARY-PART OPERATORS RE3PEC TIVELY> AND 

/ HAVE PRECEDENCE > * OR /. 

/THE FOLLOWING DECLARATIONS ARE IMPLICIT IN THE 
/COMMENTS: 

/ COMPLEX A (2**M )» U> W> T 
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0 200 
0201 
0202 
0203 
0 20 4 
0205 
0 20 6 
0 20 7 
0210 
021 1 
0212 
0213 
021 4 
021 5 
0216 
021 7 
0220 
0221 
0222 
0223 


0224 

0225 

0226 

0227 

0230 


0206 01 


0000 
0000 
40 67 
0224 
1 40 7 


/STORAGE ALLOCATION 
/ 

/AbSOLUlE SYMBOLS 
/ 



DUMMY ARGPTR /ARGUMENT POINTER 

0100 

ABSYM T 

100 

0100 

ABSYM TREAL 

100 

0101 

ABSYM T1 

101 

0101 

ABSYM AIADR 

101 

0102 

ABSYM T2 

102 

0102 

ABSYM ARVADR 

102 

0103 

ABSYM TIMAG 

103 

0103 

ABSYM IREV 

103 

0104 

ABSYM CNT 

104 

0105 

ABSYM NM2 

105 

0106 

ABSYM M 

106 /COPY OF CALLING PARAMETER 

0107 

ABSYM I SIGN 

10 7 / ” 

0110 

ABSYM N 

1 10 

011 1 

ABSYM IP 

1 1 1 

0112 

ABSYM II 

1 12 

01 1 3 

ABSYM AFIELD 

1 13 

0114 

ABSYM AADR 

1 1 4 

0115 

ABSYM L 

1 1 5 

01 1 6 

ABSYM J 

1 1 6 

0117 

ABSYM LE 

1 1 7 

0120 

ABSYM LEI 

120 


/RELOCATABLE SYMBOLS 


ARGPTR* /USED ONLY AT START OF PROGRAM 

0000 

UREAL* BLOCK 3 


0000 

0000 

0000 

UIMAG* BLOCK 3 

/REAL AND I MAG PARTS OF U 

0000 



0000 

0000 

WREAL* BLOCK 3 


0000 

0000 

0000 

WIMAG* BLOCK 3 

/REAL AND IMAG PARTS OF W 

0000 

0000 

2014 

FPONE* 2014 


0000 

0 


0000 

0 /FL PT 1 .0 


2026 

PI* 2026 


2207 

220 7 


7325 

7325 /FL PT 3 

.1415926535 

0200 01 

ADRiJ* UREAL /SINCE NON-NUMERIC LITERALS ARE 


ADRW* '/JREAL /NOT ALLOWED 
/ 

/PROGRAM STARTS HERE WITH CODE TO 
/PICK UP ARGUMENTS 
/ 

FFT* BLOCK 2 

TAD I FFT /CDF FOR A ARRAY 


01 
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0231 

31 1 3 


DCA 

AFI ELD 

0232 

2225 


INC 

FFT# 

0233 

4067 


TAD 

I FFT 

0234 

0224 

01 



0235 

1 40 7 




0236 

311 4 


DCA 

AADR /ADR OF A ARRAY 

0237 

2225 


INC 

FFT# 

0240 

4067 


TAD 

I FFT /CDF FOR M 

0241 

0224 

01 



0242 

1 40 7 




0243 

3200 


DCA 

ARGPTR 

0244 

222 5 


INC 

FFT# 

0245 

4067 


TAD 

I FFT /ADR OF M 

0246 

0224 

01 



0247 

1 40 7 




0250 

3201 


DCA 

ARGPTR# 

0251 

2225 


INC 

FFT# 

0252 

4067 


TAD 

I ARGPTR /VALUE OF M 

0253 

0200 

01 



0254 

1 40 7 




0255 

3106 


DCA 

M 

0256 

4067 


TAD 

I FFT /CDF FOR I SIGN 

0257 

0224 

01 



0260 

1 40 7 




0261 

3200 


DCA 

ARGPTR 

0262 

2225 


INC 

FFT# 

0263 

4067 


TAD 

I FFT /ADR OF ISIGN 

0264 

0224 

01 



0265 

1 407 




0266 

3201 


DCA 

ARGPTR# 

0267 

222 5 


INC 

FFT# 

0270 

4067 


TAD 

I ARGPTR /VALUE OF IS 

0271 

0200 

01 



0272 

1 40 7 




0273 

3107 


DCA 

/ 

/GET 

I SIGN 




UP DATA FIELD PART OF 


/WHICH REFER TO THE "A" ARRAY LATEH IN THE PROGRAM < 
/ 


0274 

1113 

TAD 

AFIELD 

0275 

6201 05 

DCA 

All 

0276 

3776 



027 7 

1113 

TAD 

AFIELD 

0 300 

3775 

DCA 

AR1 

0 30 1 

1113 

TAD 

AFIELD 

0302 

3774 

DCA 

AR2 

0303 

1113 

TAD 

AFIELD 

0304 

3773 

DCA 

AI2 

0305 

1113 

TAD 

AFIELD 

0 30 6 

3772 

DCA 

AIR 

0 30 7 

1113 

TAD 

AFIELD 

0310 

3771 

DCA 

AIPR 

031 1 

1113 

TAD 

AFIELD 

0312 

3770 

DCA 

All 

0313 

1113 

TAD 

AFIELD 

031 4 

3767 

DCA 

AIPI 

031 5 

1113 

TAD 

AFIELD 
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0 31 6 

3766 

OCA 

AI1R 

0 31 7 

1113 

TAD 

AFIELD 

0 320 

3765 

DCA 

AI2R 

0 321 

1113 

TAD 

AFIELD 

0 322 

3764 

DCA 

AI 1 I 

0323 

1113 

TAD 

AFIELD 

0324 

3763 

DCA 

AI2I 

0325 

1113 

TAD 

AFIELD 

0326 

37 62 

DCA 

A1 

0327 

1113 

TAD 

AFIELD 

0330 

3761 

DCA 

/ 

/ 

/ 

TAD 

A 2 

N=2**M 

0331 

1 106 

M 

0332 

4760 

JMS 

EXPON 

0333 

31 10 

DCA 

/ 

/ 

/ 

3TA 

iM 

NM2=N - 2 

0334 

7344 

CLL HAL /AC=-2 

0335 

1110 

TAD 

N 

0336 

3105 

DCA 

/ 

/ 

/ 

CLA 

NM2 

DO 7 I s 1> NM2 

0337 

7201 

I AC 

0340 

31 12 

DCA 

II 

0341 

1113 

TAD 

AFIELD /FOOL ASSEMBLER BY OVERLAYING ITS CDF 

0342 

3757 

DCA EXCHNG /WITH CDF FOR THE FIELD IN WHICH "A" RESIDED 
D07AS* /START OF DO 7 LOOP SCOPE 


/ 


/ IREV = BIT REVERSE OF I 

/ 


0343 

3103 

DCA 

IREV 

0344 

1112 

TAD 

II 

0345 

310 1 

DCA 

T1 

0346 

1 106 

TAD 

M 

0347 

7141 

CIA 

CLL 

0350 

3102 

DCA 

T2 

0351 

1101 

RVLOOP > TAD T1 

0352 

7110 

CLL 

RAR 

0353 

3101 

DCA 

T1 

0354 

1 103 

TAD 

IREV 

0355 

7004 

RAL 


0356 

5377 



0357 

0423 01 



0360 

12 60 01 



0361 

1246 01 



0362 

1236 01 



0363 

1032 01 



0364 

1026 01 



0365 

1016 01 



0366 

1012 01 



0367 

0760 01 



0370 

0 7 50 01 



0371 

0744 01 
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0 372 
0373 
0374 
0375 
0376 
0377 
0 400 
0 401 
0 402 
0 403 
0 40 4 


0 40 5 
0 406 
0 407 
0410 
041 1 


0412 
041 3 
041 4 
041 5 
041 6 
041 7 
0 420 
0 421 
0 422 
0 423 
0 424 
0 425 
0 42 6 
0 427 
0 430 
0 431 
0 432 
0 433 
0 434 
0435 


0436 
0437 
0 440 
0 441 
0 442 
0 443 
0444 
0 445 


06 
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0734 

01 



0702 

01 



0666 

01 



0646 

01 



0632 

01 



7000 




3103 


DCA 

IREV 

2102 


ISZ 

T2 

4045 


JMP 

RVLOOP 

7410 




5776 






/ 




/ 

IF (I 



/ 


1 103 


TAD 

I REV 

7041 


CIA 


1112 


TAD 

II 

7700 


SMA 

CLA 

5236 


JMP 

D0 7AF 



/ 




/ 

T = A 



/ 

A ( I ) 



/ 

ACIKE 



/ 


1112 


TAD 

II 

6201 

05 

JM3 

SUB3C 

4775 




3101 


DCA 

AI ADR 

1 103 


TAD 

IREV 

4775 


JMS 

SUBSC 

3102 


DCA 

ARVADR 

1374 


TAD 

(-6 

3104 


DCA 

CNT 

6201 

05 

EXCHNG* TAD I 

1 501 




3100 


DCA 

T 

1 502 


TAD 

I ARVADR 

3501 


DCA 

I AIADR 

1 100 


TAD 

T 

3502 


DCA 

I ARVADR 

2101 


INC 

AI ADR 

2102 


INC 

ARVADR 

2104 


ISZ 

CNT 

5223 


JMP 

EXCHNG 


2112 
1112 
7041 
1 105 
7700 
4045 
7410 
5773 


<= I) GO TO 7 


ACIREV) 
> = f 


CONTINUE 


/ 

/7 
/ 

D07AF* INC II 
TAD II 
CIA 

TAD NM2 
SMA CLA 

JMP D0 7AS /REENTER LOOP 


/ 

/NO W 


THAT BIT REVERSAL IS DONE* DIDDLE AADR 
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0 446 
0447 
0 450 


0451 

0452 


0 453 
0454 
0455 
0456 


0457 
0 460 
0461 


0462 
0463 
0 464 
0465 
0466 


0 467 
0470 
0 471 
0472 
0473 
0474 
0475 
0476 
0477 
0 500 
0 501 
0 502 
0 503 
0 50 4 
0505 
0 50 6 
0 507 
0510 
051 1 


/TO REFLECT THE FACT THAT WE WILL BE USING 
/FORTRAN-TYPE SUBSCRIPTS FROM NOW ON* WHEREIN 


1114 


/THE ARRAYS START FROM THE 
/THE ZEROTH ELEMENT. 

/ 

TAO AADR 

1 3 74 


TAD (-6 

311 4 


DCA AADR 

7201 


/ 

/ DO 20 L = 1* M 

/ 

CLA IAC 

31 1 5 


DCA L 

1115 


D020AS* /START OF SCOPE OF 
/ 

/ LE = 2**L 

/ 

TAD L 

6201 

05 

JMS EXPON 

4772 
31 1 7 


DCA LE 

1117 


/ 

/ LEI = LE / 2 

/ 

TAD LE 

7110 


CLL RAR 

3120 


DCA LEI 

7240 


/ 

/ U = <1.0* 0.0) 

/ 

STA 

1 771 


TAD ADRU 

3010 


DCA 10 

4770 


JMS ONESET 

4767 


JMS ZERSET 

1120 


/ 

/ W = COMPLX( COS(PI 

/ 

TAD LEI 

1366 


TAD (-1 

7440 


SZA 

5277 


JMP TWOTST 

4765 


JMS SET W 

4764 


JMS MONESET /IF LE1=1* W 

4767 


JMS ZERSET 

5763 


JMP WMERGE 

1366 


TWOTST* TAD (-1 

7 640 


SZA CLA 

531 5 


JMP GO TRIG 

6201 

05 

JMS SETW 

4765 

4767 


JMS ZERSET /IF LE1=2* W = 

1 107 


TAD ISIGN 

7700 


SMA CLA 

5312 


JMP PLUS 

4764 


JMS MONESET 

5763 


JMP WMERGE 


>* SINCPI/LE1) ) 


= ( -1 


0 . ) 


<0 


SGN(ISIGN>*1 *0 > 
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0512 

6201 05 

PLUS# JMS ONESET 

0513 

4770 


0514 

5763 

JMP EMERGE 

051 5 

1120 

GOTRIG# TAD LEI /CALL UPON TRIG FUNCTION 

0 51 6 

4033 

CALL 0 #FLOT 

051 7 

0002 06 


0 520 

4033 

CALL 1#STO /LEI IN FL PT 

0 521 

0103 06 


0 522 

6201 05 

ARG T 

0 523 

0100 


0 524 

4033 

CALL 1#FAD 

0525 

0104 06 


0526 

6201 05 

ARG PI 

0 52 7 

0217 01 


0 530 

4033 

CALL 1 # FO V 

0531 

0105 06 


0 532 

6201 05 

ARG T 

0533 

0100 


0 534 

4033 

CALL 1 #STO 

0535 

0103 06 


0 536 

6201 05 

ARG T /T=PI/LE1 

0537 

0100 


0 540 

4033 

CALL 1#C0S 

0541 

0106 06 


0 542 

6201 05 

ARG T 

0 543 

0100 


0544 

4033 

CALL 1#STO 

0 545 

0103 06 


0546 

6201 05 

ARG WREAL / WREAL=COS (T) 

0 547 

0206 01 


0 550 

4033 

CALL 1 * SIN 

0551 

0107 06 


0552 

6201 05 

ARG T 

0553 

0100 


0 554 

1107 

TAD ISIGN /TEST FOR SIGN REVERSAL 

0555 

7700 

SMA CLA 

0 556 

5762 

JMP WISTO 

0557 

4033 

CALL 0# CHS 

0 560 

0010 06 


0 561 

5377 


0562 

0600 01 


0563 

0604 01 


0 564 

1320 01 


0565 

1327 01 


0566 

7777 


0567 

1303 01 


0 570 

1311 01 


0571 

0222 01 


0 572 

1260 01 


0573 

0343 01 


0574 

7772 


0575 

1271 01 


0576 

0351 01 


0577 

7000 


0 600 

4033 

WISTO# CALL 1# STO 

0601 

0103 06 


0 602 

6201 05 

ARG WIMAG /WIMAG=SGNCISI GN)*SIN(T) 


T0 EVALUATE 


W 
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0 603 


0604 
0 60 5 


0 60 6 
0 607 


0610 
061 1 
0612 


061 3 
0614 
061 5 
061 6 
0617 
0 620 
0 621 
0 622 
0623 
0 624 
0625 
0626 
0627 
0 630 
0631 
0632 
0 633 
0634 
0635 
0636 
0637 
0 640 
0641 
0 642 
0643 
0644 
0645 
0646 
0647 
0650 
0651 
0652 
0653 
0654 


021 1 

01 

WMERGE* 

/ 

/ DO 20 J = 1 * LEI 

/ 

7201 


CLA IAC 

31 1 6 


DCA J 

D020BS* /START OF SCOPE OF INNER 
/ 

/ DO 10 I ■ J» N» LE 

/ 

1116 


TAD J 

3112 


DCA II 

D010AS* /START OF SCOPE OF DO 10 
/ 

/ IP = I + LEI 

/ 

1112 


TAD II 

1120 


TAD LEI 

31 1 1 


DCA IP 

/ 

/ T=A(IP) * U 

/ 

1 1 1 1 


TAD IP 

6201 

4776 

05 

JMS SUtsSC 

3247 


DCA AR1# /SET UP AKG PSEUDO OP 

1247 


TAD AK1# /ADDRESS PARTS FOR AUP 

3267 


DCA AR2# /REAL PART 

1 247 


TAD AR1# /NOW IMAG PART 

1375 


TAD C 3 

3233 


DCA All# 

1233 


TAD All# 

3303 


DCA AI2# 

4033 


CALL 0*CLEAR 

001 1 

06 


4033 


CALL 1*FAD 

0104 

06 


6211 

0000 


All* ARG 0 /-> .IMAG. A(IP) 

4033 


CALL 1 * FMP 

0112 

06 


6201 

05 

ARG UIMAG 

0203 

01 


4033 


CALL 1 * STO 

0103 

06 


6201 

0100 

05 

ARG TREAL /TR£AL=. IMAG. AUP) * i 

4033 


CALL 1 * FAD 

0104 

06 


621 l 
0000 


AR1 * ARG 0 /-> .REAL. ACIP) 

4033 


CALL 1 * FMP 

0112 

06 


6201 

05 

-ARG UREAL 

0200 

01 


4033 


CALL 1 * FSB 




PAGE 10 


-19- 


0655 

0113 

06 


0656 

6201 

05 

ARG TREAL 

0657 

0100 



0 660 

4033 


CALL 1> STO 

0661 

0103 

06 


0662 

6201 

05 

ARG TREAL / 

0663 

0100 



0664 

4033 


CALL 1> FAD 

0665 

0104 

06 


0666 

621 1 


AR2, ARG 0 / 

0667 

0000 



0 670 

4033 


CALL 1,FMP 

0 671 

0112 

06 


0672 

6201 

05 

ARG UIMAG 

0673 

0203 

01 


0674 

4033 


CALL 1>STO 

0675 

0103 

06 


0676 

6201 

05 

ARG TIMAG 

0677 

0103 



0 700 

4033 


CALL 1,FAD 

0 701 

0104 

06 


0702 

621 1 


AI2, ARG 0 / 

0 703 

0000 



0 70 4 

4033 


CALL 1,FMP 

0 70 5 

0112 

06 


0706 

6201 

05 

ARG UREAL 

0 70 7 

0200 

01 


0710 

4033 


CALL 1, FAD 

071 1 

0104 

06 


0712 

6201 

05 

ARG TIMAG 

0713 

0103 



0714 

4033 


CALL 1,STO 

0715 

0103 

06 


0716 

6201 

05 

ARG TIMAG / 

0717 

0103 


/ 

/ A (IP 

/ 

0 720 

1247 


TAD AR1 # /- 

0721 

3345 


DCA AIPR# 

0 722 

1233 


TAD All# 

0 723 

3361 


DCA AIPI# / 

0724 

1112 


TAD II 

0725 

4776 


JMS 5UBSC 

0 726 

3335 


DCA AIR# /- 

0727 

1335 


TAD AIR# 

0 730 

1375 


TAD <3 

0731 

3351 


DCA All# 

0732 

4033 


CALL 1,FAD 

0 733 

0104 

06 


0734 

621 1 


AIR, ARG 0 7 

0735 

0000 



0736 

4033 


CALL 1,FSB 

0737 

0113 

06 


0 740 

6201 

05 

ARG TREAL 

0741 

0100 



0 742 

4033 


CALL 1,STO 


A < IP ) = A(I) - T 
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0 743 
0744 
0745 
0746 
0 7 47 
0750 
0751 
0752 
0753 
0754 
0755 
0756 
0757 
0 760 
0761 


0762 

0763 

0764 

0765 

0766 

0767 

0770 

0771 

0772 

0773 

0774 

0775 

0776 

0777 

1000 

1001 

1002 

100 3 
1004 
1 00 5 
1006 
1 007 
1010 

101 1 
1012 

1013 

1014 

1015 
101 6 
1017 
1020 
1 021 
1022 

1023 

1024 
1 025 
1026 
1027 
1 030 


1 1 


0103 

06 



621 1 
0000 


AIPR, 

ARG 0 

4033 

0104 

06 

CALL 

1 * FAD 

621 1 
0000 


All* 

ARG 0 

4033 

0113 

06 

CALL 

1 * FS8 

6201 

0103 

05 

ARG 

TIMAG 

4033 

0103 

06 

CALL 

1 * STO 

621 1 
0000 


AIP I * 

ARG 0 


/ 

/ 

/ 


A (I ) = AID + T 


1335 


TAD AIR# 

6201 

3774 

05 

DCA AI IK# 

1 33 5 


TAD AIR# 

3773 


DCA AI2R# 

1 351 


TAD All# 

37 72 
5377 


DCA AI1I# 

1027 

01 


1017 

01 


1013 

0003 

01 


1271 

7000 

01 


4045 

7410 

1 776 


TAD All# 

3233 


DCA AI2I# 

4033 


CALL 1 * FAD 

0104 

06 


6201 

0100 

05 

ARG TREAL 

4033 


CALL 1 * FAD 

0104 

06 


621 1 
0000 


AI1R, ARG 0 

4033 


CALL 1 * STO 

0103 

06 


621 1 
0000 


AI2R, ARG 0 

4033 


CALL 1 * FAD 

0104 

06 


6201 

0103 

05 

ARG TIMAG 

4033 


CALL 1 * FAD 

0104 

06 


621 1 
0000 


AI11, ARG 0 

4033 


CALL 1,STO 
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1 031 

0103 06 


1032 

621 1 

AI2I> ARG 0 /-> .IMAG. Ad) 

1 033 

0000 

/ 

/10 CONTINUE 

/ 

001 OAF * /END OF DO 10 LOOP SCOP; 

1 034 

1112 

TAD II 

1035 

1117 

TAD LE 

1036 

3112 

DCA II 

1037 

1112 

TAD II 

1040 

7041 

CIA 

1041 

1 110 

TAD N 

1042 

7700 

SNA CLA 

1043 

4045 

JMP D010AS /LOOP AGAIN 

1044 

7410 


1045 

5775 

/ 

/ U = U * W 

/ 

1046 

4033 

CALL 1 * FAD 

1047 

0104 06 


1050 

6201 05 

ARG UIMAG 

1 051 

0203 01 


1052 

4033 

CALL 1 *FMP 

1053 

0112 06 


1054 

6201 05 

ARG WIMAG 

1055 

0211 01 


1056 

4033 

CALL 1»STO 

1057 

0103 06 

- 

1060 

6201 05 

ARG TREAL /TREAL=UIMAG* WIMAG 

1061 

0100 


1062 

4033 

CALL 1>FAD 

1063 

0104 06 


1064 

6201 05 

ARG UREAL 

1065 

0200 01 


1066 

4033 

CALL 1»FMP 

1067 

0112 06 


10 70 

6201 05 

ARG WIMAG 

1071 

0211 01 


1072 

4033 

CALL 1»STO 

1073 

0103 06 


1074 

6201 05 

ARG TIMAG /TIMAG=UREAL* WIMAG 

1 075 

0103 


1076 

4033 

CALL 1»FAD 

1077 

0104 06 


1 100 . 

6201 05 

ARG UREAL 

1 101 

0200 01 


1 102 

4033 

CALL 1»FMP 

1 103 

0112 06 


1 104 

6201 05 

ARG WREAL 

1 105 

0206 01 


1 106 

4033 

CALL 1»FSB 

1 107 

0113 06 


1 1 10 

6201 0 5 

ARG TREAL 

1111 

0100 


1112 

4033 

CALL 1»STO 
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1113 

0103 

06 


1114 

6201 

05 

ARG UREAL /UREAL=UREAL* WREAL - UIMAG* WIMAG 

1115 

0200 

01 


1116 

4033 


CALL 1»FAD 

1117 

0104 

06 


1 120 

6201 

05 

ARG UIMAG 

1121 

0203 

01 


1 122 

4033 


CALL 1 * FMP 

1 123 

01 12 

06 


1 124 

6201 

05 

ARG WREAL 

1 125 

0206 

01 


1 126 

4033 


CALL 1»FAD 

1 127 

0104 

06 


1 130 

6201 

05 

ARG TIMAG 

1 131 

0103 



1 132 

4033 


CALL 1*ST0 

1 133 

0103 

06 


1 134 

6201 

05 

ARG UIMAG /UIMAG=UREAL* WIMAG + UIMAG* WREAL 

1 135 

0203 

01 



/ 

/20 CONTINUE 

/ 

D020BF, /END OF SCOPE OF INNER DO 20 LOOP 


1136 

211 6 

INC 

J 


1 137 

1116 

TAD 

J 


1 1 40 

7041 

CIA 



1 1 41 

1120 

TAD 

LEI 


1 142 

7700 

SMA 

CLA 


1 1 43 

4045 

JMP 

D02QBS 

/LOOP AGAIN 

1 144 

7410 




1 145 

5774 






D020AF » /END OF SCOPE OF 

1 1 46 

21 1 5 

INC 

L 


1 147 

1115 

TAD 

L 


1 150 

7041 

CIA 



1 151 

1 106 

TAD 

M 


1 152 

7700 

SMA 

CLA 


1 1 53 

4045 

JMP 

D020AS 

/LOOP AGAIN 

1 1 54 

7410 




1 1 55 

5773 

/ 





/ 

IF 

CISIGN >= 0 ) 



/ 



1 1 56 

1107 

TAD 

ISI GN 


1 1 57 

7710 

SPA 

CLA 


1 1 60 

5363 

JMP 

NORMLZ 


1 161 

4040 

RETRN FFT 


1162 

0001 06 

/ 





/ 

M s 

2*N 



/ 



1 1 63 

1110 

NORMLZ# TAD 

N 

1 164 

7104 

CLL 

RAL 


1 1 65 

3106 

DCA 

M 




/ 





/ 

T = 

1 . / N 
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1 1 66 

1 1 10 


1 1 67 

4033 


1 1 70 

0002 

06 

1 171 

5377 


1 1 73 

0453 

01 

1 174 

0606 

01 

1 175 

0610 

01 

1 1 76 

0751 

01 

1 1 77 

7000 


1 200 

4033 


1 201 

0103 

06 

1 202 

6201 

05 

1 203 

0103 


1 20 4 

4033 


1205 

0104 

06 

1206 

6201 

05 

1 207 

021 4 

01 

1 210 

4033 


121 1 

0105 

06 

1212 

6201 

05 

1 213 

0103 


1214 

4033 


1215 

0103 

06 

1216 

6201 

05 

121 7 

0100 



TAD N 

CALL 0»FLOT 


CALL 1 * STO 
ARG TIMAG 
CALL 1> FAD 
ARG FPONE 
CALL 1,FDV 
ARG TIMAG 
CALL 1 > STO 
ARG T 


1220 7201 

1221 3112 


1222 1114 

1223 1377 

1224 3114 


1225 

1112 


1 226 

7104 


1 227 

1112 


1 230 

1114 


1231 

3237 


1 232 

1237 


1 233 

3247 


1 234 

4033 


1 235 

0104 

06 

1 236 

621 1 


1237 

0000 


1 240 

4033 


1 241 

01 12 

06 


/ 

/ DO 40 I = 1# M 

/ 


CLA IAC 
DCA II 

/NOW WE WILL SEQUENCE THROUGH THE A 
/IF IT WERE A REAL ARRAY INSTEAD OF 
/THUS WE NEED TO DIDDLE THE ADDR OF 
TAD AADR 
TAD (3 
DCA AADR 

DO 40 AS > /START OF SCOPE OF DO 40 LOOP 
/ 

/ AAU) * AA( I ) * T 

/ 

/WHERE 
/ 

/ 

/ 


ARRAY JUST AS 
COMPLEX. 

A AGAIN 


THE FOLLOWING IS 
EQUIVALENCE CAC1 
REAL AA 


ASSUMED 
)> AA <1 ) ) 


TAD II 

CLL RAL /COMPUTE SUBSCRIPT 
TAD II 
TAD AADR 
DCA Al# 

TAD Al# 

DCA A2# 

CALL 1»FAD 


Al, ARG 0 /-> AA<I) 
CALL 1»FMP 
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1 242 
1 243 
1244 
1 245 
1 246 
1247 


1 2 50 
1 251 
1 252 
1253 
1 254 
1255 
1 256 
1 257 


6201 0 5 
0100 
4033 
0103 06 
621 1 
0000 


2112 
1112 
7041 
1106 
7700 
5225 
40 40 
0001 06 


ARG T 

CALL 1»STO 
A2> ARG 0 /-> AA(I) 

/ 

/ 40 CONTINUE 

/ 

D040AF, /END OF SCOPE OF DO 40 LOOP 
INC II 
TAD II 
CIA 
TAD M 
SMA CLA 

JMP DO40AS /LOOP AGAIN 
RETRN FFT 


/ 

/END OF COMPUTATIONAL PART OF SUBROUTINE 
/UTILITY PROGRAMS FOLLOW 
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/UTILITY SUBROUTINES FOR FFT SUBROUTINE 
/ 



/RAISE 2 TO THE 

INTEGER POWER IN THE AC 



/RETURN INTEGER 
/ 

ANSWER IN AC 

1260 

0000 

/ 

EXPON* 0 


1 261 

7041 

CIA 


1262 

3100 

DCA T 


1 263 

7001 

I AC 


1 264 

7104 

EXL* CLL RAL 


1265 

2100 

ISZ T 


1 266 

5264 

JMP EXL 


1 267 

•6201 05 

JMP .1 EXP ON 


1 270 

5660 




/ 


/COMPUTE THE ADDRESS IN THE A ARRAY 
/CORRESPONDING TO THE INTEGER SUBSCRIPT 
/IN THE AC. RETURN ADDRESS IN AC 
/ 


1271 

0000 

SUBSC* 0 

1272 

3302 

DCA 

SUTEMP 

1 273 

1302 

TAD 

SUTEMP 

1 274 

7104 

CLL 

RAL 

1 275 

1302 

TAD 

SUTEMP 

1 276 

7104 

CLL 

RAL 

1 277 

1114 

TAD 

A A DR 

1 300 

6201 05 

JMP 

I SUBSC 

1 301 

5671 



1 302 

0000 

SUTEMP* 0 


/ 


/MOVE FL PT ZERO INTO WRDS -> TO BY AXIO 
/ 


1 303 

0000 

ZERSET* 

0 

1 304 

6201 05 

DCA I 

10 

1 305 

3410 



1306 

3410 

DCA I 

10 

1 30 7 

3410 

DCA I 

10 

1310 

5703 

JMP I 

ZERSET 


/ 

/MOVE FL PT 1.0 INTO WRDS -> TO BY AXIO 
/ 


1 31 1 

0000 

ONESET* 0 

1 312 

6201 0 5 

TAD FPONE 

1 313 

1 776 


1 31 4 

3410 

DCA I 10 

1 31 5 

3410 

DCA I 10 

1316 

3410 

DCA I 10 

1 317 

57 1 1 

JMP I ONESET 


/ 


/MOVE FL PT -1.0 INTO WRDS -> TO BY AXIO 
/ 


1 320 

0000 

MONESET 

* 0 

1 321 

1375 

TAD 

(601 4 

1 322 

6201 0 5 

DCA 

I 

10 

1 323 

3410 




1 324 

3410 

DCA 

I 

10 

1 325 

3410 

DCA 

I 

10 
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1 326 

5720 


JMP I 
/ 

MONESET 




/SET 

/ 

AX REG 10 

1 327 

0000 


SETW, 

0 

1 330 

7240 


STA 


1 331 

6201 

05 

TAD 

ADR W 

1 332 

1 774 




1 333 

3010 


DCA 

10 

1 334 

5727 


JMP 

I SET W 

1 374 

0223 

01 



1 375 

601 4 




1 376 

0214 

01 



1 377 

0003 





TO POINT TO W 


END 




